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I.  INTRODUCTION 


Vector  transformations  between  the  mean-of-date  and  instantaneous  earth-fixed  reference 
systems  are  performed  through  application  of  the  g  and  £  transformation  matrices,  or  their  trans¬ 
poses,  to  the  vector  of  interest  (these  matrices  are  described  in  detail  in  the  next  section).  These 
transformation  matrices  are  functionally  dependent  upon  the  nutation  in  longitude,  A\p;  the  nuta¬ 
tion  in  obliquity,  Ae;  and  the  equation  of  the  equinoxes,  AH.  The  A Ac,  and  AH  angles  are  time 
dependent  and  can  be  computed  in  a  straightforward  manner  at  each  time  for  which  they  are  re¬ 
quired.  However,  since  each  individual  computation  is  quite  lengthy,  such  an  approach  is  im¬ 
practical  if  g  and  £  matrices  are  required  for  a  great  many  times. 

A  more  practical  approach  can  be  used  to  generate  the  required  A\p,  Ae,  and  AH  angles  for 
situations  in  which  many  g  and  £  matrices  are  needed.  Required  values  for  these  angles  are  ob¬ 
tainable  by  interpolating  between  values  which  have  been  precomputed  at  discreet  times  over  the 
timespan  of  interest.  However,  the  interpolation  introduces  errors  into  A^,  Ae,  and  AH,  which  are 
manifested  as  errors  in  the  associated  transformed  vector. 

The  following  sections  of  this  work  discuss  in  detail  these  manifested  errors  in  the  transformed 
vector,  as  well  as  the  magnitude  of  the  interpolation  errors  in  A Ae,  and  AH.  In  the  next  section  a 
general  first-order  error  equation  is  developed  that  analytically  describes  the  upper  error  bound  on 
the  vector  moduli  obtained  from  mean-of-date/instantaneous  earth-fixed  vector  transformations  due 
to  small  errors  in  the  A^,  Ae,  and  AH  angles.  This  analytical  result  is  applied  to  a  six-point  Lagrange 
interpolator  in  section  III  where  numerical  data  are  used  to  verify  the  accuracy  of  the  analytical 
result,  as  well  as  assess  the  computational  quality  of  the  interpolator  for  this  application.  The  last 
section  summarizes  the  results  and  conclusions  drawn  from  this  analysis. 


II.  ERROR  EQUATION 


A  first-order  analytic  expression  for  the  upper  bound  of  the  error  introduced  into  a  vector 
modulus  by  the  coordinate  transformation  from  the  mean-of-date  reference  system  to  the  instan¬ 
taneous  earth-fixed  system  is  developed  in  this  section.  In  this  discussion  it  is  assumed  that  this 
error  is  an  interpolation  error  induced  by  using  transformation  matrix  elements  computed  from 
angular  arguments  obtained  by  an  interpolation  process. 

The  transformation  of  a  vector  in  the  mean-of-date  reference  system,  ^ ,  to  a  vector  in  the 
true  of  date  system,  Xp  is  accomplished  by  the  following  application  of  the  nutation  matrix,  £: 

•  (1) 

Similarly,  the  transformation  of  to  a  vector  in  the  instantaneous  earth-fixed  system,  £g,  is  per¬ 
formed  through  the  rotation, 


1 


t 


*E  *  fi*r  ■ 


(2) 


where  g  is  the  earth-fixed  transformation  matrix.  Using  Equation  (1)  in  Equation  (2)  provides  the 
complete  mean  of  date  to  instantaneous  earth-fixed  transformation  given  by 


It  should  be  mentioned  that  both  g  and  £  are  orthogonal  transformations,  i.e. 
g'1  *  g* 

and 

£‘  *  £' 
so  that 

=  C  *T  • 

*r  =  • 

and 

-  G‘  B  • 

In  the  above  expressions  the  superscript  “t”  means  matrix  transpose. 

The  g  and  £  matrices  have  the  following  forms: 

cos  A  sin  A  0 
g  *  -sin  A  cos  A  0 

0  0  1 


(3) 

(4) 

(5) 

(6) 

(7) 

(8) 


(9) 


and 

[cosA^  -sin  A^  cos  7  -sin  sin  ? 

cos  e  sin  A^  cosecosAiJ’Cosf+sincsin?  cosecosA^sinl-sinecos?  J,  (10) 
sin  e  sin  A^  sin  e  cos  A^  cos  7  -  cos  e  sin  7  sin  e  cos  A^  sin  7  +  cos  e  cos  7 

where  A  is  the  longitude  of  the  Greenwich  meridian  from  the  true  vernal  equinox  at  time  t  given  by 
A  -  H,,*  AH(t)  +  w(t-Atj)  ;  (1!) 
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V 


A^  is  the  nutation  in  longitude;  Ae  is  the  nutation  in  obliquity;  ?  is  the  mean  obliquity  of  date;  and 
e  is  true  obliquity  of  date  given  by 

e  =  e  +  Ae  (12) 

The  terms  in  Equation  (11)  have  the  following  definitions: 

Hq  =  mean  hour  angle  of  Greenwich  at  zero  hours  universal  time  of  the  day  of 
interest 

AH(t)  =  equation  of  the  equinoxes  =  Atf/  cos  e 

g>  *  mean  sidereal  rotation  rate  of  the  earth 

At2  =  correction  for  irregular  earth  rotation 

Suppose  that  the  angles  A Ae,  and  AH  are  obtained  from  interpolations  between  A Ac,  and 
AH  values,  which  have  been  precomputed  at  discreet  times  over  the  timespan  of  interest,  and  that 
these  interpolations  introduce  the  errors  fi^,  fie,  and  6H  into  their  values,  i.e., 


A A^0  +  S\l>  ,  (13) 

Ae  =  Ae0  +  fie  ,  (14) 

e  =  ?  +  Ae0  +  fie  *  eQ  +  fie  ,  (IS) 

and 

A  -  A0  +  fiH  ,  06) 

where  A^0,  Ae«,  e„,  and  AQ  are  the  “em>rless”values  of  A^,  Ae,  e,  and  A,  respectively.  Substi¬ 
tuting  the  last  four  equations  into  the  expressions  for  g  and  assuming  that  the  errors  are  very 
small;  and  using  the  small  angle  approximations 

cos  fid  *  1,  sin  fid  *  fid  for  fid  «  1  ,  (17) 

provide  the  following  first-order  results: 

ft  »  Ro+  «Hft,  ,  (18) 

and 

£  -Co+  8*£t+  *£2  >  <19> 


"herefto  and  Co  are  the  “errorless”  transformation  matrices  given  by  Equations  (9)  and  (10)  with 
A,  A^,  and  e  replaced  with A *0,  and  e0,  respectively.  The  error  matrices  have  die  following 
forms: 
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-  sin  A, 


ft,  =  -  cos  A0  -  sin  A0 


sin  A^0 


cos A^0  cos? 


cos  A^0  sin  ? 


£,  =  I  cose0cosA^0  -sin A^0  cose0  cos?  -sin  A\P0  cos e0  sine 


sin  e0  cos  A^0  -  sin  A^0  sin  e0  cos  ?  -  sin  A^0  sin  eQ  sin  e 


Q2  =  -sine0sinA^0 -sin  e0  cos  A0O  cos?  +  cos  e0  sin? -sin  ^  cos  Arsine -cos  cos?  .(22) 
cosCgSinA^Q  cos  eQ  cos  A^0  cos  ?  +  sin  eQ  sin  ?  cos  cos  Arsine- sin  eQ  cos? 


Equations  (18)  and  (19)  may  be  used  in  Equation  (3)  to  give  the  following  first-order  result: 


*E  =  (  Bo  Co  +  8*  Bo  Cl  +  *€  Q2  +  SH  Co  )  c 

Assume  that  Xc  can  be  decomposed  into  the  sum  of  an  “errorless”  part,  XF  ,  and  a  part  corre- 
*■  -*  o 

sponding  to  the  interpolation  induced  error,  SX£.  Since 

'  xe0=  feo  Go  ’  * 


«*e  “  (  8*  fto  G,  +  «e  Ro  G2  +  «H  fti  Co) 


Taking  the  norm  of  both  sides  of  Equation  (25)  gives 


I8*eI  -  I  (  8*  ft0  £,  +  5e  Ro  Q2  +  «H  ftj  Go) 


Upon  application  of  the  Schwartz  and  triangle  inequalities,  the  last  equation  may  be  rewritten  as 
the  following  inequality: 


-jjpj-  <  WllftoCil  +  IfiellftoGjI  +  ISHMft,  Col 


The  matrix  norms  |  ft,  £j  |  have  the  values  given  by1  : 


1  Noble,  B., 


,  Inc.,  Englewood  Cttfft,  New  Jeney,  1969)  pp.  427-431. 


I 


w 


ifcCjf  “  *S'L  ,  (»- (01,02.10) 


where 


\imx  *  max2mum  eigenvalue  of  (J^  Gj)11  (ft,  £j)  .  (2S 

The  superscript  “h”  in  the  last  equation  means  hermitian  transpose.  Since  R,and  Gj  are  real 
matrices,  the  hermitian  transpose  is  simply  the  matrix  transpose,  so  that 

(S*  G/  (Si  Gj)  =  %  Gj)‘  (Bi  Cj) 

(3( 

-  cia/BiCi  • 

For  the  case  i  =  0,  Equation  (4)  applies,  and  Equation  (30)  becomes 

(Ro  GjV  <fto  Gj)  •  G|  q  0  -  1,2).  (31 

The  maximum  eigenvalue  X^.  is  then  obtained  from  the  solution  of  the  characteristic  equation 

d®t  { Gj  Gj  -  *ojI}  =  o,o  =  1,2)  ,  (3: 

where  “det  (^ }”  means  the  determinant  of  and X is  the  identity  matrix.  Using  Equations  (21) 
and  (22)  to  form  £*  Gj  0  =  1  >  2)  gives 


-2  sin  &i)0  cos  Atf>0  cos  e  -2  sin  A^„  cos  A^0  sin  e 


Gj  G1  *  "2  sin  A^0  cos  Atlig  cos  ? 

-2  sin  A^0  cos  A^0  sin  t 


cos2  6 


sin  ?  cos  f 


sin  e  cos  e 


sin2  A^„  sin  A^„  cos  A^0  cos  I  sin  Ajl»0  cos  A^0  sin  ? 

Gj  Cj  *  *in  A#0  cos  A^0  cos  ?  cos2  A^0  cos2  ?  +  tin2  ?  -tin2  Atf^  tin  ?  cos  I  .(j 

sin  A#0  cot  A^0  sin  ?  -tin2  A^0  sin  ?  cos  ?  cos2  A^0  tin2  ?  +  cos2  ? 

Upon  substitution  of  tbs  taut  two  expressions  into  Equation  (32)  and  after  much  algebraic 
manipulation,  it  is  found  that  the  associated  characteristic  equations  reduce  to  the  simple  forms: 

*oi  [  *oi  -  2*oi  "  <4sin*  A*0coo*A*0  -  1)]=  0  (3 


S 


and 


C 


i 


I 


! 


\2  [xoa  -  2xo2  +  l]  “  0  .  (36) 

which  yield 

Xo,  =  {0,  1  ±2sinAtf0  cosA0o}  ,  (37) 

and 

*02  =  {0,  1,  1}  .  (38) 

Thus 


*01  max  =  1  +  2  I  sin  A*0  cosA*0  !  ,  (39) 

and 

*02  max  *  1  •  (40) 

Using  Equation  (30)  it  is  seen  that  for  i  =  1,  j  =  0,  the  maximum  eigenvalue  XJ0  is  ob¬ 
tained  from  the  solution  of  the  characteristic  equation  m&' 

det  {G*  ft,  Go  ~  *,oX}  *  0  •  (4D 

However,  by  inspection,  it  is  found  that 

ft‘,  ft,  =  X  -  X  >  (42) 

where 

(-i)fck  =  ^k3  •  (43) 

In  the  last  expression,  6mn  is  the  Kronecker  delta.  Using  this  result,  the  properties  of  the  identity 
matrix  and  the  relation  (see  Equation  (5)) 

Go  Co  =  X  (44) 

allow  Equation  (41)  to  be  rewritten  as 

det(GS[(l>XIO)X-X]Coj  *  0  •  (45) 
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Application  of  the  properties  of  determinants  to  this  expression  allows  it  to  be  recast  into  the  form 


det  {Go}  det  {(1  -  X, 0)X  -  i  }  det  {  Go  }  =  0  .  (46) 

Since 

**lci I- 1«  {«■)-  tojcT)  •  (47) 

then  Equation  (46)  reduces  to 

det  {(1  -X10)X-i.  }=  0  .  (48) 

The  solution  to  this  equation  is  straightforward,  and  it  is  found  that  the  characteristic  equation 
reduces  to 

X,0  (Xio  -  1)(X10  -  D  «  0  I  (49) 

yielding 

Xl0  =  {0,  1,  1  }  .  (50) 

Thus 


^1 0  max  *  1  •  (5D 

Use  of  Equations  (28),  (39),  (40),  and  (51 )  in  Equation  (27)  gives  the  final  form  of  the  error 
equation: 

I  a*E  i 

- -  <  (1  +  2  1  sin  A*0  cos  AiK  l)*  I  5^  |  +  |  5e  I  +  1  8H  1  .  (52) 

l  x»  I 

Consider  now  the  error  associated  with  the  inverse  transformation,  i.e.,  that  given  by  Equation 
(8).  Inserting  Equations  ( 1 8)  and  (19)  into  Equation  (8),  retaining  only  first-order  terms,  using  the 
inverse  of  Equation  (24),  and  applying  the  Schwartz  and  triangle  inequalities  to  the  resulting  norms 
provide  the  inequality 

I«*mI 

—  <  I  **  I  I  Q\  fto  I  +  I  6e  i  1  Gl2  fto  I  +  I  5H  I  I  GUI  I  •  (53) 

I  XE  | 

Analagous  to  Equation  (28),  the  matrix  norms  I  G{  B]  I  have  the  values  given  by: 


where 


-fy  =  maximum  eigenvalue  of  (g*  g1,)h  (gj  g*) , 


(55) 


i.e.,  the  is  the  maximum  eigenvalue  obtained  from  the  solution  of  the  characteristic 

J,max 

equation 

det  "  >iiX)  =  0  .  (56) 

For  the  case  when  (j i)  =  (01 )  the  transpose  of  Equation  (44)  applies  and  Equation  (56)  reduces  to 


det  {  Bi  B}  "  70,  i  }  =  0  • 

However,  it  is  seen  from  Equation  (20)  that 

B,  B‘,  =  Bi  B,  . 

so  that  Equation  (42)  applies  and  Equation  (57)  reduces  to  Equation  (48).  Thus 
7oi  ”7jo  ~  1  • 

max  max 

For  the  cases  when  (ji)  =  ( 1 0,  20)  use  may  be  made  of  the  fact  that 

Bo  Bo  =  X 

to  rewrite  Equation  (56)  as 

det|Bo[£j£j  -  7j0x]gl)J=  0  •  *  i>2) 


or 


^  {Bo}  det  -  7j0 X)  det  { ft*  }  =  0,0=  1,2) 


Since 


then  Equation  (62)  reduces  to 

<fctjc,q-Ti0i|-o. 


(57) 


(58) 


(59) 


(60) 


(61) 


(62) 


(63) 


(64) 
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Using  Equation  (2 1 )  it  is  found  that  Qj  Qj  is  given  by  Equation  (33)  with  e  replaced  everywhere 
with  e0  and  Equation  (64)  reduces  to  an  equation  of  the  same  form  as  Equation  (35)  since  it  is  in¬ 
dependent  of  I  and  e0.  Therefore 

7l0m«  =  X01ro«  =  1  +  2  I  sin  A*o  008  A^o  *  '  (65) 

Similarly,  using  Equation  (22),  it  is  found  that 

£2  £2  =  i  ”  &  >  (66) 

where 

=  5JU  5ki  •  (67) 


As  before,  is  the  Kronecker  delta.  Since  this  matrix  is  diagonal,  the  eigenvalues  lie  along  the 
diagonal  and  are  given  by  Equation  (38)  so  that 


^20. 


■  \ 


02. 


=  1 


(68) 


Therefore  the  right-hand  side  of  inequality  (53)  is  the  same  as  that  of  inequality  (52),  i.e.,  the  upper 
error  bounds  are  identical  for  the  transformation  and  its  inverse,  and  inequality  (52)  is  valid  as 
written  or  with  the  ‘E’  and  ‘M’  subscripts  interchanged. 


in.  NUMERICAL  APPLICATION:  THE  SIX-POINT  LAGRANGE  INTERPOLATOR 


In  order  to  examine  the  utility  of  the  error  equation  that  is  derived  in  the  previous  section,  a 
numerical  test  case  is  discussed  in  this  section.  Since  the  six-point  Lagrange  interpolator  is  fre¬ 
quently  used  as  the  method  for  computing  interpolated  nutation  angles,  it  has  been  chosen  for  use 
in  this  application.  This  method  involves  fitting  a  polynomial  of  degree  five  to  six  of  the  data 
points  of  the  type  being  fitted  (i.e.  A<£,  Ac,  AH)  and  centering  the  fit  such  that  three  of  the  points 
lie  on  either  side  of  the  interpolation  point.  The  spacing  between  all  data  points,  i.e.,  the  inter¬ 
polation  interval,  is  assumed  to  be  equal  and  uniform.  The  polynomial  f(t)  used  to  interpolate  the 
nutation  angle  data  has  the  following  form: 


t2  <  t  <  t 


3  ’ 


where 


and  f(t)  represents  Atfr(t),  Ae(t),  or  AH(t). 


(69) 


(70) 
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The  primary  objectives  to  be  satisfied  by  the  analysis  of  this  numerical  application  are 

(i)  Numerically  determine  the  variation  of  the  maximum  interpolation  errors  I  6  $  1 _ . 

•  5e  I  max* and  •  SH  Lax  with  interpolation  interval  for  the  six-point  Lagrange  interpolator 

(ii)  Numerically  validate  inequality  (52) 

(iii)  Quantify  the  maximum  errors  induced  by  this  interpolation  method  over  a  range  of 
Earth  satellite  position  vector  moduli 

To  generate  data  for  this  analysis,  values  for  the  errors  B\p(t),  8e(t),  and  6H(t)  resulting  from 
the  six-point  Lagrange  interpolation  algorithm  were  computed  every  15  min  using  a  1-day  interpo¬ 
lation  interval  over  the  30-day  span  from  5  January  1984  through  3  February  1984.  The  variations 
of  A\£/,  Ae,  and  AH  during  this  time  are  shown  in  Figure  1.  The  maximum  errors  |  Bip  |_r .  |  fie  1^, 
|  5H  lmax  for  this  timespan  were  isolated  and  subjected  to  recomputation  using  0.75-,  0.50-,  and 
0.25-day  interpolation  intervals.  These  error  computations  were  performed  using  sections  of  A.  R. 
Darnell’s  ROTATE  computer  program2  (i.e.,  subroutines  SMTRAN,  NOD,  and  COMPHO)  to  gener¬ 
ate  the  “errorless”  values  A\p0,  Ae0,  and  AHg  at  the  discreet  time  points  required  by  the  interpolator, 
as  well  as  the  values  of  these  quantities  at  the  interpolation  times  t.  The  “errorless”  and  interpo¬ 
lated  Aij/,  Ae,  and  AH  values  were  subtracted  to  obtain  the  fi^,  fie,  and  fiH  errors,  respectively,  at 
each  interpolation  time.  It  should  be  noted  that  all  computations  performed  by  the  ROTATE 
program  are  done  in  a  fashion  which  is  consistent  with  the  new  J2000  fundamental  reference  sys¬ 
tem  and  associated  transformations  that  are  planned  for  use  in  1984.  Figure  2  shows  the  resulting 
variation  of  |  fi^»  lmax,  |  fie  lmax,  md  |  fiH  imax  with  interpolation  interval. 

The  ROTATE  program  was  also  used  to  evaluate  Equation  (3)  for  an  arbitrary  of 
( 1 000, 1 000, 1 000)  over  an  1 8-day  span  from  1 7  January  1 984  through  3  February  1 984.  The  g 
and  Q  matrices  were  evaluated  using  both  “errorless”  and  interpolated Jup,  Ae,  and  AH  values,  and 
the  respective  XE  vectors  were  subtracted  to  obtain  the  error  vector  5XE  induced  by  the  inter¬ 
polation  errors  B\p,  fie,  and  8H.  This  process  was  repeated  every  1 5  min  using  a  1-day  interpo¬ 
lation  interval.  The  maximum  daily  errors  |  5  \p  lDtyi  ,  I  fie  lD,yj _ ,  I  6H  lDtyi _ ,  and 

^  '  max  *  max  nux 

•  sxe  b*yimgx  0  =  1 7>  18,...,  34)  were  isolated  and  used  to  validate  inequality  (52)  written 
in  the  form 

I  Ww  5  j  I  H  lD.„_  -  I  *  *  I «» lD.,w  | 

where  use  has  been  made  of  the  fact  that  Av('0  is  small  so  that 

(1  +  2  |  sin  A^0  cos  A^0  |)w  *  1  .  (72) 

These  results  are  shown  in  Figure  3.  There  the  dots  represent  the  common  logarithm  of 
the  daily  maximum  error  of  I  fij£E  |  computed,  using  the  ROTATE  program;  the  X*s  represent 
the  common  logarithm  of  the  right-hand  side  of  inequality  (71),  using  the  daily  maximum  values 
for  |  B\p  |,  |  Be  I,  and  |  fiH  |  computed  by  ROTATE;  and 

1^1  =  1  (1000, 1000,  1000)  I  «  1732.05  .  (73) 


I  X*,  I  .  (71) 


2Dimefl,  A.  R.,  Tmafor  nation  Between  the  ECI  w d  ECEF  CoortjMW  SytfiM.  NSWC  TN  (Dakfem,  VhgMa,  IN 2) 
In  ptepariUon. 
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LEGEND 

X  -  MAXIMUM  DAILY  ERROR 
BOUND  FROM  EQ.  (71) 

•  -  MAXIMUM  DAILY  ERROR 
COMPUTED  USING  ROTATE 


X  «  *  ,  «  « 


•  •  ■ 


DAY/1984 


Figure  3.  Maximum  Daily  Error  Induced  in  X£  by  Using  Interpolated 
Nutation  Angles  in  the  ft  Q  Transformation  Matrices 
for  |  |  =  1732.05 


As  can  be  seen  from  this  figure,  the  computed  error  is  bounded  in  every  case  by  the  right-hand 
side  of  inequality  (71).  Thus  inequality  (71)  is  valid  over  this  timespan. 

To  satisfy  objective  (iii),  the  following  assumption  is  made:  the  maximum  interpolation 
errors  I  6^  lmiX,  I  de  lm#x,  and  |  SH  |mgx  obtained  during  the  5  January  1984  through  3  February 
1984  timespan  are  typical  of  those  obtained  from  a  six-point  Lagrange  interpolator  for  any  time- 
span  when  using  that  nutation  model,  which  is  consistent  with  the  J2000  fundamental  reference 
system  and  associated  transformations.  From  this  assumption,  inequality  (52)  may  be  rewritten 
to  represent  the  maximum  error  bound  on  I  |  induced  by  the  use  of  nutation  angles  computed 
from  a  six-point  Lagrange  interpolator.  This  inequality  is  given  by 


s|i 


l*«l. 


where  use  has  been  made  of  approximation  (72).  The  results  for  the  maximum  interpolation 
errors  shown  on  Figure  2  can  be  used  to  evaluate  the  factor  in  curly  braces  in  the  last  expression 
for  any  interpolation  interval  between  1 .00  and  0.25  days.  The  result  obtained  for  a  1-day  inter¬ 
polation  interval  when  1 5£M  |  ■  1732.05  is  shown  by  the  hatched  line  on  Figure  3.  Note  that 
this  result  bounds  all  the  daily  computed  errors  and  daily  error  bounds  computed  from  in¬ 
equality  (71). 


f 
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